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CN ; Abstract 

Freezing of simple liquid metals and the relative stabilities of competing crys- 
^ . talline solids are investigated using thermodynamic perturbation theory, the inter- 

actions between ions being modeled by effective pair potentials derived from pseu- 
dopotential theory. The ionic free energy of the solid phase is calculated, to first 
order in the perturbation potential, using classical density-functional theory and an 
accurate approximation to the hard-sphere radial distribution function. Free energy 
calculations for Na, Mg, and Al yield well-defined freezing transitions and structural 
free energy differences for bcc, fee, and hep crystals in qualitative agreement with 
'^ \ experiment. 
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1 Introduction 



A microscopic description of freezing/melting phenomena presents a formidable 

challenge to present-day condensed matter theory [1]. Among the simple liq- 
Q\ ' uids, liquid metals have received limited attention in this field, stemming from 

^ ' the relative complexity of their interatomic potentials and thermodynamic 

'j^ . properties [2]. Previous related studies include those of Stroud and Ashcroft 

[3], who addressed the melting of simple metals using variational principles; 
,.^ , Igloi et al. [4] , who studied the freezing of liquid aluminum and magnesium by 

qj I means of density-functional (DP) methods; and Moriarty and coworkers [5], 

O • who calculated the structural stabilities of third-period simple metals using 

first-principles techniques. 

/\ ■ In this contribution, we report on a new study of the freezing phenomena of 

c^ ■ simple metals that exploits recent advances in thermodynamic perturbation 

theory for solids. In line with Ref. [4], and most modern theories of freezing, 
our approach incorporates classical DP theory [6] in calculating the Helmholtz 
free energy. Although formally rigorous, DP theory is limited in practice by 
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insufficient knowledge of the free energy functional for most except hard-body 
systems. To circumvent this limitation, we apply Weeks-Chandler-Andersen 
(WCA) perturbation theory [2,7], which splits the effective interatomic pair 
potential into reference and perturbative parts and maps the reference system 
onto an effective hard-sphere (HS) system. For the inhomogeneous solid phase, 
DF theory provides an adequate approximation for the HS free energy, while 
an accurate model of the HS radial distribution function yields directly the 
first-order perturbation free energy. The homogeneous liquid phase is treated 
within the same perturbative framework, using essentially exact expressions 
for the HS free energy and pair distribution function. This treatment of the 
liquid and solid by the same theoretical approach constitutes an important 
requirement for a consistent description of phase coexistence. 

The next section outlines the theoretical methods used to calculate effective 
pair potentials and free energies of simple metals. In Sec. 3 we present and 
discuss predictions of the theory for freezing transitions and relative stabilities 
of competing crystal structures. Finally, in Sec. 4 we summarize and conclude. 



2 Theoretical Methods 



Restricting attention to the simple metals, we construct effective density (p)- 
dependent interatomic pair potentials, (f){r;p), and a volume-dependent, but 
structure insensitive, contribution to the binding energy, Ev{p), via pseud- 
potential theory, assuming linear response of the conduction electrons to the 
ions. In practice we use empty-core pseudopotentials and the Ichimaru-Utsumi 
electron dielectric function, which are known to give reliable results for bulk 
liquid properties [8]. Taking 0(r; p) and Ev{p) as the basic microscopic input, 
we proceed to calculate the thermodynamics of the liquid and solid phases 
by means of WCA perturbation theory. Thus, we separate 0(r; p) at its ffist 
minimum into a short-range, purely repulsive, reference part, 0o(^); and a 
long-range, weakly oscillating, perturbative part, 0p(r;p), and map the refer- 
ence system onto an effective HS system with a HS diameter prescribed by 
the well-known WCA criterion. This separation allows the ionic contribution 
to the Helmholtz free energy to be calculated within the coupling-constant 
formalism [2]. To ffist order in 0p(r; p), the resulting general (inhomogeneous) 
expression for the total Helmholtz free energy reads 

F[p{r)] = Fns[p{r)] + -^ J dr' r" gus{r'; [p(r)]) 0p(r';p) + Ey{p), (1) 





where N is the number of particles, V the volume, and Fhs[p(i")] and gns{f^] [p{^)]) 
the free energy and radial distribution function (RDF), respectively, of the HS 



reference system, both functionals of the equihbrium one-particle number den- 
sity p(r). The RDF of the sohd is defined as an orientational and translational 
average of the two-particle density [9]. For the liquid we use the homogeneous 
version of Eq. (1), taking for the HS free energy and RDF the Carnahan- 
Starling and Verlet-Weis parametrizations of computer simulation data [2] . 

To calculate the free energy of the HS solid we invoke classical DF theory 
[6]. The DF approach is based on the existence of a functional jF[p(r)] of the 
density p(r) that satisfies a variational principle, according to which jF[p(r)] 
is minimized - for a given average density and external potential - by the 
equilibrium density, its minimum value equaling the Helmholtz free energy 
F. In practice, jF[p(r)] is split into an (exactly-known) ideal-gas contribution 
jFid[p(r)] and an excess contribution jFex[p(r)], the latter depending entirely 
upon internal interactions. Here we approximate jFcx[p(r)] by the modified 
weighted-density approximation (MWDA) [10]. This approximation maps the 
excess free energy per particle of the solid onto that of a corresponding uniform 
fluid of effective density p, according to 

•^e™^[p(r)] = iV/Hs(p), (2) 



where the effective density, defined as 

P = ^y dry dr'p(r)p(r')w;(|r-r'|;p), (3) 

is a self-consistently determined weighted average of p(r). The weight function, 
w(r), is specified by normalization and by the requirement that ^,!^^°^[p(r)] 
generates the exact two-particle direct correlation function in the uniform limit 
(see Ref. [10] for details). 

Practical calculation of jFHs[p(r)] and 5'hs(^; [p(i^)]) requires specifying the solid 
density, ie., the coordinates of the lattice sites and the shape of the density dis- 
tribution about these sites. Here we consider the hep, fee, and hee crystals with 
the density distribution modeled by the usual Gaussian ansatz, introducing a 
parameter a determining the width of the distribution. This parametrization 
allows the ideal contribution to the free energy functional to be accurately 
approximated by J-'id/N = (3/2)A;BTln(Q;A^) —5/2, where A is the thermal de 
Broglie wavelength, and yields the HS free energy as the minimum with respect 
to a of the approximate functional jFHs[p(r)] = jFjd[p(r)] + jF^^°'^[p(r)]. 

For the RDF appearing in the perturbation contribution to the free energy 
[Eq. (1)] we adopt the approach of Rascon et al. [9]. Expressing 5'hs(^; [p(i')]) 
as a sum over coordination shells, this scheme approximates the distributions 
of the second and higher coordination shells in simple mean-field fashion, but 
parametrizes the first-shell distribution in such a way as to incorporate nearest- 



neighbour correlations. The first-peak parameters are determined by imposing 
sum rules, e.g., the virial theorem and coordination number, so as to specify 
the contact (r = d) value and shape (area and first moment) of the first peak 
(for details see Ref. [9]). For a given solid structure, the lattice distances, the 
coordination number, the a that minimizes jFjjg[p(r)], and the HS pressure, 
derived directly from FHs[p(r)], then combine to determine 5'hs(^5 [p(i^)]) ^'^d 
hence the perturbation free energy. 

Following the above procedure, we obtain the ionic free energy, -Fion, which is 
the only stracture-dependent contribution to the total free energy, and which 
suffices for assessing relative stabilities of different solid structures. To fully 
describe freezing/melting, however, it is essential to consider the total free 
energy, F = Fion + Ey, where Ey is an additional contribution from the con- 
duction electrons. The electronic free energy (or volume energy), Ey, although 
independent of structure, depends on the average density and thus influences 
the densities of coexisting phases. Here we take for Ey the usual expression 
following from second-order perturbation theory [8] and consistent with the 
effective pair potential. 



3 Results and Discussion 



Figure 1 illustrates the density dependence of the effective HS diameter, d, 
for three different crystal structures of Al near its melting point. According 
to the WCA prescription, d depends both on the reference pair potential and 
on the flrst peak of the HS pair distribution function. Since the flrst peak 
is identical for hep and fee crystals, d is the same for these two close-packed 
structures. Differences in the flrst peak for the bee crystal result in a smaller 
diameter for that more open structure. Because of the sensitivity of Fhs to the 
HS packing fraction, these structural variations in d, though only of the order 
of 1 %, can amount to signiflcant variations in the free energy. In passing, 
we note that the widely-used Barker-Henderson prescription for d [2] neglects 
such structural dependencies and that only knowledge of the HS RDF permits 
consistent application of the more accurate WCA prescription [9]. 

Predicted structural free energy differences for Na, Mg, and Al are displayed in 
Fig. 2. Evidently the theory correctly predicts the observed equilibrium struc- 
tures for these three metals, namely bcc-Na, hcp-Mg, and fee-Al, although the 
density at crossover from fee- to bee-Na is somewhat overestimated. Compared 
with predictions of flrst-principles approaches, our results agree qualitatively 
well. For Mg and Al the stability trends of Fig. 2 match those from general- 
ized pseudopotential theory (GPT) and the linear muffin-tin orbitals (LMTO) 
method [5]. As well, for Na we predict the same trend as GPT, while LMTO 
fails for this case. Compared with available experimental data, our structural 



free energy differences are in somewhat better quantitative agreement than 
those of either GPT or LMTO. For example, for Mg and Al we predict dif- 
ferences between hep and fee free energies of —1.6 and 3.2 mRy, respectively, 
compared with —0.6 and 1.7 mRy from GPT, —0.7 and 1.6 mRy from LMTO, 
and —1.5 and 4.2 mRy from experiment (extrapolations of thermochemical al- 
loy data, as quoted in Ref. [5]). It should be emphasized, however, that while 
the first-principles methods calculate ground-state (T = 0) energies, thermo- 
dynamic perturbation theory yields free energies of finite-temperature sys- 
tems. We note further that application of perturbation theory is restricted to 
those structures for which the HS system is at least metastable, i.e., for which 
jFHs[p(r)] has a local minimum. This excludes, e.g., the diamond structure of 
Si, the HS diamond crystal being unstable. 

Finally, we investigate the freezing/melting transition by comparing predicted 
free energies for the liquid and solid phases. Figure 3 exemplifies the com- 
parison for Al near its observed melting temperature. Crossing of the ionic 
free energy curves for the liquid and fee phases implies a distinct transition 
near the observed equilibrium solid density. Similar comparisons for Na and 
Mg also demonstrate liquid-solid phase transitions in qualitative agreement 
with observation. In this respect, the effective-pair-potential approach reason- 
ably describes the many-body interacting system. Note that the liquid-solid 
crossover density serves as a useful upper (lower) bound on the liquid (solid) 
density at coexistence. Adding the volume-energy contribution - identical for 
the two phases - while not affecting the crossover density, does alter the shape 
of the curves. As a result, the densities of coexisting liquid and solid phases, 
as determined by a Maxwell common tangent construction, depend on Ey. 
An accurate volume energy is known to be essential for calculating realistic 
pressures and bulk moduli of metals from pseudopotential theory at the level 
of effective pair interactions [11,12]. Similarly, we find the liquid-solid coex- 
istence analysis to be sensitive to the form of Ey, particularly its curvature 
with respect to p. In fact, the sensitivity is such as to preclude here a reliable 
determination of coexistence densities. Given the reasonable predictions for 
ionic free energies, this apparent limitation of the theory may reflect a need to 
go beyond the linear response approximation to include higher-order response 
effects in the volume energy (see also [8], p. 42ff). The relative significance 
of such effects, as well as distinctions [12] between the real-space representa- 
tion of Ey used here and the reciprocal-space representation of Ref. [3], are 
interesting issues deserving of further attention. 



4 Conclusions 



Summarizing, working within an effective-pair-potential framework, we have 
implemented a consistent form of thermodynamic perturbation theory, one 



that is of comparable accuracy for both the hquid and sohd phases, and used 
it to study structural stabilities and freezing behaviour of third-period sim- 
ple metals. Predictions for structural free energy differences are at least as 
accurate as those from two first-principles methods and are in generally good 
qualitative agreement with experiment. On the basis of ionic free energies, 
freezing/melting transitions are predicted near the observed densities. As the 
densities of coexisting liquid and solid phases are found to depend sensitively 
on the electronic free energy, a complete coexistence analysis has not been at- 
tempted pending more accurate knowledge of this part of the free energy. The 
computationally practical approach demonstrated here should prove useful in 
future studies of the freezing of liquid metal mixtures, as well as in guiding 
more sophisticated ab initio simulations. 
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Fig. 1. WCA effective liard-sphere diameter in atomic units (Bolir radii) for Al at 
temperature T = 930 K in hep, fee, and hce crystal structures. 



Fig. 2. Free energy differences between various crystal structures vs. reduced atomic 
volume, where f^o is the observed equilibrium atomic volume at atmospheric pres- 
sure: (a) Na at T = 300 K; (b) Al and Mg at T = 930 K. 



Fig. 3. Ionic free energy vs. average density for liquid and /cc-solid phases of Al at 
T = 930 K. 
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Fig. 2a: Denton, Kahl and Hafner 
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Fig. 2b: Denton, Kahl and Hafner 
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Fig. 3: Denton, Kahl and Hafner 
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